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Abstract 

We study the half-plane problem for the elastic wave equation subject to a 
free surface boundary condition, with particular emphasis on almost incompressible 
materials. A normal mode analysis is developed to estimate the solution in terms of 
the boundary data, showing that the problem is boundary stable. The dependence 
on the material properties, which is difficult to analyze by the energy method, is 
made transparent by our estimates. The normal mode technique is used to analyze 
the influence of truncation errors in a finite difference approximation. Our analysis 
explains why the number of grid points per wave length must be increased when 
the shear modulus (/i) becomes small, that is, for almost incompressible materials. 
To obtain a fixed error in the phase velocity of Rayleigh surface waves as — > 0, 
our analysis predicts that the grid size must be proportional to /i 1 / 2 for a second 
order method. For a fourth order method, the grid size can be proportional to /i 1//4 . 
Numerical experiments confirm these scalings and illustrate the superior efficiency 
of the fourth order method. 



1 Introduction 

Consider the half-plane problem for the two-dimensional elastic wave equation in a ho- 
mogeneous isotropic material. By scaling time to give unit density, the displacement with 
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Cartesian components (u, v) T is governed by 

u tt = (iAu + (A + n)(u x + v y ) x + Fi(x,y,t), 

x > 0, — oo < y < oo, t > 0, (1) 

v tt = [iAv + (A + ii)(u x + v y )y + F 2 (x,y,t), 

where (F±, F 2 ) T is the internal forcing. Here, A and \i > are the first and second Lame 
parameters of the material. We assume that both parameters are constant and A > 0. 
The displacement is subject to initial conditions 

u(x,y,0) = f w (x,y), jv(x, y, 0) = / u (a;, y), 

\ x > 0, —00 < y < 00. (2) 

u t (x,y,0) = f 2 o{x,y), [v t (x,y,0) = f 2 i{x,y), 

In this paper we consider normal stress boundary conditions along the x = boundary, 
u x + i 2 v y = gi(y,t), 

x = 0, —00 < y < 00, t > 0, (3) 

u y + v x = g2\y,t), 

where g\ and g 2 are boundary forcing functions, and 

A 

7 - 



2/i + A 



When gi = and g 2 = 0, ([2]) is called a free surface boundary condition. 
Since time was scaled to give unit density, the elastic energy is given by 



00 ^00 



E{t) = 1-1 I (u 2 t + vf) + X(u x + Vyf + n (2u 2 x + 2v 2 + (u y + v x f) dx dy. (4) 

2 J -oo Jo 

It is well known (see e.g. Achenbach pQ, pp. 59-61) that the elastic energy satisfies 

^ /»oo /*oo 

= J J {u t F 1 + v t F 2 ) dxdy 

POO 

(u t ((2/i + \)u x + Xvy) + vtn(u y + v x ))\ x=Q dy. 

In particular, without boundary and interior forcing, the elastic energy is conserved, 

E(t) = E(0), t>0, g 1 = g 2 = 0, F 1 = F 2 = 0. (5) 

Note that the elastic energy is a semi-norm of the solution. The energy estimate bounds 
this semi-norm in terms of the initial data and the internal forcing (F 1? F 2 ) T . For this rea- 
son, the elastic wave equation is a well-posed problem. However, the energy estimate does 
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Figure 1: A contour plot of a Rayleigh surface wave as function of (x,y) at t = for 
a material with A = 1 and = 0.01. The w-component is shown to the left and the v- 
component to the right. The contour levels are given between -0.5 and 0.5, with spacing 
0.05. Red and blue lines correspond to negative and positive values, respectively. The 
zero level is plotted in black. 

not provide detailed insight into how the solution depends on the material parameters, or 
the boundary data. 

The material parameters, in particular the ratio /i/A, strongly influences the accuracy 
of numerical solutions of the elastic wave equation. As a motivating example, we propagate 
a Rayleigh surface wave using a second order accurate finite difference method. In the 
numerical experiment, we make the y-direction 1-periodic and take the wave length to be 
one. A free surface boundary condition is imposed at x — 0. The Rayleigh surface wave 
propagates harmonically in the y-direction and decays exponentially in x, see Figured! We 
take A = 1 and vary /i, which gives the surface wave a phase velocity that is proportional 
to JJl. We discretize the elastic wave equation on a grid with grid size h, corresponding 
to P = 1/h grid points per wave length; further details of this numerical experiment are 
presented in § [51 In Figure [21 we report the error in the numerical solution at time t = 20. 
For the smaller values of /i, a large number of grid points per wave length are needed 
to obtain an acceptable error level and a second order convergence rate. For the finest 
mesh with 200 grid points per wave length, the error increases by more than an order of 
magnitude (from 3.76 • 10 -3 to 5.09 • 10~ 2 ), when \i decreases by two orders of magnitude 
(from 10 _1 to 10~ 3 ). Note that the gradient of the exact solution only depends weakly 
on fi and is of the order 0(1/2tt) for all values of /i > 0. Hence, the loss of accuracy is 
not due to poor resolution in space. Furthermore, the phase velocity of the surface wave 
becomes slower and slower as /i — > 0, while the time step is governed by y/X + 3/z, which 
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Figure 2: Max error, normalized by the max norm of the exact solution, at time t = 20 
in the numerical solution of the Rayleigh surface wave problem when A = 1 and different 
values of ft. The error is shown as function of the number of grid points per wave length, 
P = l/h. 



tends to vA = 1. Hence, the temporal resolution of the surface wave only improves as 
fi -»• 0. 

In this paper we use a normal mode analysis to explain the loss of accuracy as fi — > 0, 
which corresponds to the incompressible limit of an elastic material. The normal mode 
analysis allows us to estimate the solution in terms of the boundary data, and makes the 
dependence on the material parameters transparent. We show that the solution is strongly 
boundary stable, except in the vicinity of the generalized eigenvalues corresponding to 
surface waves. Here the solution is as smooth as the boundary data, i.e., only boundary 
stable (see [5] for definitions of these stability concepts). We develop a modified equation 
model of the truncation errors in the numerical calculation, where we view the discretized 
boundary conditions as a perturbation of the exact boundary conditions. This analysis 
reveals how perturbations of the boundary conditions influence the solution, and how the 
material parameters enter in the relation. 

To analyze the solution of (CO)-®, we follow the technique used by Kreiss, Ortiz and 
Petersson [5] and split the problem into two parts. First we consider a Cauchy problem, 
where the definition of the forcing and the initial data are extended to the whole of M 2 (x). 
Secondly, we subtract this solution from the solution of the half-plane problem to obtain 
a new half-plane problem, where only the boundary data do not vanish. This is a very 
natural procedure because all the difficulties and many physical phenomena arise at the 
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boundary. The new half-plane problem is analyzed in detail using the Fourier-Laplace 
transform method, leading to estimates of the solution in terms of the boundary data. 

The remainder of the paper is organized in the following way. The properties of 
the Cauchy problem are briefly discussed in Section [2j The normal mode analysis of 
the half-plane problem is developed in Section [3] We discuss the eigenvalue problem 
in Section I3.1ti3.2[ leading to necessary conditions for a well-posed problem. Boundary 
estimates are derived in Section 13.31 In Section HI we use the normal mode theory to 
perform a modified equation analysis of the discretized boundary conditions. This analysis 
shows how the number of grid points per wave length must be increased to maintain a 
given error level in the numerical solution when // — > 0. For the second order method, 
the grid size must be proportional to /i 1//2 , while it suffices to take h ~ /i 1//4 for the fourth 
order method. These scalings are confirmed by the numerical experiments in Section [5j 
illustrating that the fourth order method is significantly more efficient than the second 
order approach, in particular for small values of /i. Conclusions are given in Section 



2 The Cauchy problem 

In this section we consider the Cauchy problem for ([I])-©. The definitions of the forcing 
functions and the initial data can be smoothly extended to the whole of M 2 (a;). For 
simplicity we use the same symbols for the extended functions as for the original ones. 

We start by deriving an equation for the divergence of the displacement, 5 = u x + v y , 
by forming the divergence of ([T]). This gives 

S u = (A + 2/x)A5 + G(x,y,t), -oo < (x, y) < oo, t > 0, (6) 

where the forcing is G = dF\/dx + dF 2 /dy. The divergence, 5 = 8(x,y,t), is subject to 
initial conditions 

o{x,y,0) = + d t (x,y,0) = + , -oo < (x, y) < oo. (7) 

By first solving the wave equation for the divergence, we can (in principle) treat the 
divergence as a forcing in the Cauchy problems for u and v, 

u tt = fiAu + F 1 ( y x,y,t), 

- oo < (x,y) < oo, t > 0, (8) 

v tt = fxAv + F 2 (x,y,t), 
where 

Fx(x,y,t) = (A + /i)5 x + i r i(x,?/,t), F 2 (x,y,t) = (A + fi)6 y + F 2 {x,y,t). 

Since 5, u, and v all satisfy scalar wave equations, we conclude that the Cauchy problem 
for the elastic wave equation is well-posed. 
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Note that the wave propagation speed in the wave equation for the divergence is 
y/X + 2/i. For G = 0, admits plane wave solutions of the type 



S(x,y,t) = e 



: iu){xznc p t) 



\J\ + 2/i. 



Hence, a wave with angular frequency £ = uc p has wave length 




(9) 



Note that L p stays bounded for A = const., /i — > 0. By taking the curl of (DO), we can also 
derive a scalar wave equation for the curl of the displacement, where the wave propagation 
speed is y/Ji. Hence, the elastic wave equation also admits plane waves with wave length 



The length of these waves tend to zero as /i — > 0. 

3 The half-plane problem 

We are interested in solutions with bounded L 2 -norm and therefore we assume 



Throughout the remainder of the paper, s = rj + z£ denotes a complex number where 
7], £ are real numbers. As a preliminary, we define the branch cut of y/a + ib by 



where a and b are real numbers, 

3.1 A necessary condition for well-posedness, the eigenvalue 
problem 

We start with a test to find a necessary condition such that the half-plane problem is well 
posed. 

Lemma 1. Let F\ = = and g% — g-i — 0. The problem (JT]) -(J3j) is not well-posed if 
we can find a non-trivial simple wave solution of the type 




for every fixed t. 



(10) 



7r < arg (a + ib) < ir, arg V a + ib 




u = U(x)e st+iulv , v = V(x)e st+iuy , 
\U\oo < oo, (V^loo < oo, Re(s) > 0, uj real. 



(11) 
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Proof. If we have found such a solution, then 

ui = U(ax)e sat+iujay , vi = V(ax)e sat+iway , 

is also a solution for any a > 0. Since Re(s) > 0, we can find solutions that grow 
arbitrarily fast in time. The problem is therefore not well-posed. □ 

We shall now discuss whether there are such solutions. Introducing ( II ip into ([T]) gives 

(s 2 + fiu 2 )U - (2fi + \)U XX - i(X + im)ojV x = 0, ^ > q 
(s 2 + (2/i + A)w 2 ) V - /iV^ - z(A + fi)uU x = 0, 

To derive boundary conditions for U and V, we insert f fTTj) into ([3]), 

x = 0. (13) 

iw£/ + V x = 0, 

Equation (fl2j) is a system of linear ordinary differential equations with constant coef- 
ficients. It can be solved using the ansatz 

U(x) = u e~ KX , V{x) = v e~ KX . (14) 

Inserting (JT3J) into ffl2l) gives a linear system for (tio,t>o) T , which can be written 

(s 2 + /iw 2 — (2/i + A)k 2 ) m + i(\ + fi)uK Vq = 0, (15) 
i(X + /j)ojku + (s 2 + (2/i + X)uj 2 - fin 2 ) v = 0. (16) 

Let 

/■ 2 , 2 2 
(, = S + fiUJ — [IK . 

Then we can write (IT5"1)-([T6"1) as 



(C- (X + fi)K 2 )u + i(X + fi)coKv = 0, (17) 
+ fi)ojKu + (C + (A + /i)w 2 )w = 0. (18) 

This system has a non-trivial solution if and only if its determinant is zero, 

[C - (A + /i)/t 2 ] [C + (A + /iV 2 ] + (A + /i) 2 u;V = 0. 

There are two possibilities. Either ( = 0, or ( + (A + [i){uj 2 — k 2 ) = 0, corresponding to 



S 2 / S 2 

k = ±\Iuj 2 -\ , or, K = ±\ UJ 2 + 



fi' V (2/i + A) 



In appendix |A] we shall prove that there is a constant 5 > such that 



Ro ( \h' 2 + — J >SRe(s), Re ^y^ 2 + A + 2jU I - rt ~ R °(- s )- IM*) > 0. 
Thus, for Re(s) > 0, there are two solutions that have bounded L 2 -norm: 

^ ) =e-^u 1 + e-^u 2 , Ui = ( U ° J | ., = 1.2. (19) 
V(x)J \v oj 

with 



«, = ,/<"■ + £. K2 = y^ + _i!_ (20) 

It is convenient to calculate the eigenvectors by inserting (ki,Ui) into (ITS1) and (k 2 ,u 2 ) 
into (fI7j). 

z(A + jj)u)KiUoi + (A + /i)w 2 f i = 0, 

-(A + [l)u) 2 U Q2 + i(X + /j)uJK 2 V 02 = 0- 

Therefore, 

«01 = W01, V 02 = «02- 

OJ K 2 

We summarize these results in the following lemma. 

Lemma 2. Assume Re(s) > and uj ^ 0. T/ien /?i 7^ k 2 and the general solution of the 
ordinary differential equation [IE) can be written as 

^ j = W01 ( + u 02 [ L I e"^, (21) 

where K\ and k 2 are en 6u (1201) . 

Remark 1. Inserting [21\) into < f77]j shows that all simple wave solutions satisfy 

is 2 

v x -u y = —u 01 e st+ ^-^ x , (22) 

«* + «y = ~ , x «02 e st+ ^-^. (23) 

(A + 2/j,)k, 2 

Hence, uqi and uq 2 are proportional to the curl and divergence of the solution, respectively. 
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Introducing (12TI) into the boundary conditions (|T3|) gives 

(1 - 7 2 )/«i/t 2 u m + (k 2 . - 7 V) n 02 = 0, (24) 
{k\ + u 2 )u m + 2tu 2 u 02 = 0. (25) 

The linear system fl24^) - (!25|) has a non-trivial solution if and only if its determinant is 
zero, 

A =: 2u 2 (l - 7 2 )K!fi: 2 - («| - 7 V) (k 2 + w 2 ) = 0. (26) 
Since 1 — 7 2 = 2/i/(A + 2/i), we can write ( 1261) in the form 

where 

^==^^ + ^-(1 + 1)', 1=^. (27) 

Note that the zeros of the determinant ( J26l) are the solutions of <p(s) = 0. 

Lemma 3. Assume w^O. TTie function </?(s) does not /iai>e any zeros /or Re (s) > 0. 

Proof. Assume there was a solution of (p(s) = with Re (s) > 0. It would correspond 
to a non-trivial solution (1401,^02) of ( 1241) - (1251) . There would therefore be a simple wave 
solution ( ITT]) where £7(x) and V(x) are given by (T2TT) . This simple wave solution would 
have Re(s) = |a;| v //IRe(s) > 0, and for this reason, its elastic energy P| would grow 
exponentially in time. However, this is contradicted by the energy estimate which says 
that the elastic energy must be constant in time. There can therefore be no simple wave 
solutions for Re(s) > 0, and the function (p(s) can not have any zeros for Re (s) > 0. □ 

As a consequence of this lemma, 

Theorem 1. The elastic wave equation (T2])-([2P with F\ = F 2 = and g\ = g 2 = 0, has 

no simple wave solutions of the type other than the trivial solution u = v = 0. 

Because ffT2]) - ffT3]) define an eigenvalue problem, we can also phrase the theorem as 
Theorem 2. The eigenvalue problem (1121) - 0131) has no eigenvalues with Re(s) > 0. 

3.2 Generalized eigenvalues 

We shall now calculate the generalized eigenvalues, i.e., roots of the determinant (I2T1) in 
the limit Re (s) — > 0+. We need to discuss 5 = 0;, £ real, and the zeros are given by 

^-^•^-^-(l-f) 2 ^, ? = (28) 

We have 
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Lemma 4. Equation l[28\) has the solution £ = ; and exactly two solutions s~o — =t^£o 
with < £o < 1- There are no solutions with £ 2 > 1. 

Proof. Inserting £ = into fl28l) shows that it is a solution. Clearly, there are no solutions 
for 1 < £ 2 < (2/i + A)//i because the first square root is purely imaginary and the 
second square root is real. Also, the second term in ip is always real and negative. For 
£ 2 > (2/i + A)//i, both square roots are purely imaginary and their product is real and 
negative. Hence both terms in <p are real and negative. We conclude that there are no 
solutions for £ 2 > 1. 

To analyze < £ 2 < 1, we denote A = X/fi and observe that the function 



has the same roots as <p. It has the properties 
1. 



#)) = 0, v(i) = ~<o, 

lb 



2. 



that is, 



dip /da =: = - (l + J + + 2 f 1 - J 

V 2 + A/ 2+A V 2 



// n 1 + A . . , 2 + 3A 

^ = — ^ > 0, ^1 = < o. 

2 + A 8 + 4A 



that is, 



^"(0) = -i±^ < 0, ^» = for - = 1 ± „ 

1 ; 2 + A 1 ^ 2 3(2 + A) 

Thus ip"(cr) has at most one sign change in < a < 1. Properties 1-3 show that ■0' has 
one sign change and the lemma follows. □ 

In Table [1] we have calculated the scaled generalized eigenvalues Sq = — £o f° r some 
values of A//a. Note that all values remain bounded in the limit A/// — > oo, i.e. /i — > 
when A = const. 

Differentiating ( 1271) gives 
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A//i 


s 2 - -£ 2 
^0 ~~ so 


kio/M 









-0.7639 


0.4858 


0.7861 


0.6036 


1 


-0.8452 


0.3933 


0.8474 


1.0610 


4 


-0.8877 


0.3350 


0.9230 


1.6045 


8 


-0.8991 


0.3175 


0.9539 


1.8360 


oo 


-0.9126 


0.2955 


1 


2.1936 



Table 1: Coefficients in the solution at the generalized eigenvalues s = ±i-JJl\cj\^ , for 
some values of A//x. 

Because <p(so) = 0, ( 127|) gives 



2/i + A s 



<~2^ 



2/i + A 2/x + A 



-2 1 + 



S2\ 3 



; — Co; 



(30) 



where C = C (A///) is real. Since s is purely imaginary, <^'(so) is also purely imaginary. 
We report numerical values of |y/(s )| in Table HJ demonstrating that <£>'(sq) is bounded 
away from zero for all values of X/fi > 0. Therefore, (p(s) has a first order zero at the 
generalized eigenvalues s"o = ±i£o- 

To calculate the eigenfunctions corresponding to the generalized eigenvalues so = ±«£o> 
we consider the two boundary conditions ( |T3|) . Evaluating the general solution (j2Tj) gives 



icuU + V x = icu ( ( 1 -\ — \ ) uqi + 2u, 



'02 



At the generalized eigenvalues, 



x = 0. 



K? = W 2 (l + S 2 ,) 



Hence, iu;[/ + 14 = if 

(2 - g)«oi + 2m 02 = 0. (31) 

If relation ( I3T1) is satisfied, also U x + ry 2 u;V = 0. The eigenfunction corresponding to 
so = ±i£o is therefore given by 



1 

^10 



_l_ i(£ 2 _ 2)e =t ^ * +iajy_K2oa; 




(32) 
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where 

I i A £2 I i /i £ £o 

These eigenfunctions, also known as Rayleigh waves (see e.g. Achenbach pQ, §5.11), rep- 
resent surface waves that propagate in the positive or negative y-direction. 

Now we consider the potential generalized eigenvalue s = 0. Relations ( 12?)) and 
(1291) show that both </j = and dip/ds = for s = 0. Differentiating ( 129]) shows that 
d 2 ip/ds 2 7^ for s = 0. Thus </?(s) has a zero of order two at s = 0. However, for s — >■ 0, 
(1201) show that both K\ — > \u\ and k% — > \u\. In this limit, boundary conditions (I24p and 
(1251) give 

«oi + u 02 = 0, s -)■ 0. 

Expanding the general solution (|2ip around 5 = shows that the eigenfunction vanishes 
identically in this limit. Thus s = is not a generalized eigenvalue. 



3.3 Boundary forcing 

As we discussed in the introduction, we split the solution of the half-plane problem ([I])-© 
into a Cauchy problem and a new half-plane problem, where only the boundary data do 
not vanish. Hence, the Cauchy problem satisfies the initial conditions and the interior 
forcing function. Its solution drives the solution of the new half-plane problem through 
a modified boundary forcing function. For example, when the half-plane problem (jTJ)-(J3J) 
has an interior forcing function with compact support in Cl, the solution of the Cauchy 
problem consists of waves propagating outwards from fl. The gradient of these waves 
along x = enter in the boundary forcing functions for the new half-plane problem. 

The estimates obtained in this and the following sections are expressed in Fourier- 
Laplace transformed space. It is clear that all these estimates have their counterpart in 
physical space. To understand the relation between both types of estimates, we refer to 
chapter 7.4 of [3 J or chapter 10 of [2]. 

We consider CQ)-© with homogeneous initial data and internal forcing, F\ = F% = 0. 
We Laplace transform the problem with respect to t, Fourier transform it with respect to 
y, and denote the dual variables by s and u>, respectively. Here u is a real number and s 
is complex. We obtain, 



(s 2 + fiu 2 )u — (2/i + \)u xx — + fj)v x = 0, 
(s 2 + (2/i + \)u 2 )v — fiv xx — + y)u x = 0, 

subject to the boundary condition 



x > 0, Re(s) > 0, (33) 



u x + i-f 2 uv = gi(u,s), 
. , . . . , v x = 0. (34) 

iuu + v x = g2{u, s), 
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Note that (it, v) T satisfy the same differential equation as (U, V) T in ( fl2|) . By LemmaEJ 
the general solution is of the form (12~T]) . i.e., 



u(x) 
v(x) 



«0! iK , e~™ + u 02 ltu (35) 



In the following, we assume u ^ 0. The case — > will be studied separately in 
appendix [Bl 

By inserting (1351) into boundary condition (134|) . we get 

(l - 7 2 ) K X K 2 Uoi + («| - 7 2 ^ 2 ) M 2 = ~«2 9l, 

[k\ + w 2 ) M i + 2u 2 uq 2 = —iu fa- 

This system corresponds to (r24l) - (T25"l) with an inhomogeneous right hand side. In terms 
of the scaled variable s defined by (T2T1) . 



fiS 2 



«! = \cu\VT+P, k 2 = \u\\\1 + ^-^-. (36) 



After some algebra, the system for (i/oi, Wo 2 ) T becomes 



(37) 



1 + «oi+«o 2 = - — • (3? 



The determinant of (I37j) - (l38l) is 



~9 / ~2\ 2 



where the function ip(s) was previously defined by (12T1) . To solve the system, we eliminate 
m 2 from (138]) and insert in ()3j 



/~x . (A + 2^i / s 2 /x z£ 2 / S 2V , 

U01 = — V 1 + 2^Ta + 2^( V 1 + Y ) - (39) 



Inserting this expression into (1381) gives 



(A + 2/x)^i / ^s 2 / § 2 \ i<? 2 /—^ L , 
M ° 2 = 2^| V *+2? V + " ^"^V 1 + XT^- (40) 
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Hence, the system (137|) - (!38|) becomes singular exactly at the roots of (p(s) = 0. For 
Re(s) > 0, Lemmas |3] and H] prove that this can only happen at the generalized eigenvalues. 
The general theory of [5] tells us that, away from the generalized eigenvalues, |</9(s)| _1 is 
bounded and the problem is therefore strongly boundary stable. 

We want to estimate the solution on the boundary in terms of the boundary forcing. 
For x = 0, the general solution (|35|) satisfies 



6(0) 
0(0) 



u 01 + u 02 , 



-Uqi 



«"2 



-U 2 



VTT" s 2 u Q i 



lid 

\cj\ 



1 



IJLS 2 

X + 2/j, 



-1/2 



(41) 



U02, 



We now discuss how the solution behaves close to the generalized eigenvalues So = 
±i£o- By Lemma HI we have < £o < 1 and both K\ and k 2 are real. Since p(So) = 0, 
Taylor expansion gives 



<p(s) = {S- so)v'(so) + 0{\s - s | 2 ). 



(42) 



Formula (J3QD and Table Q] shows that \(p'(s )\ > C ~ 0.6 for all X/fi > 0. We have 
s — sq = rj + i(£ — £o), an d to leading order in < rj <C 1, 



ip{S) 



So 



which leads to the estimates 



s |>?7, \<p{s)\> 



y/fi\u 



■Co, v > 0. 



(43) 



(44) 



For rj > 0, the system (|371)-(l38l) is non-singular and we can substitute ( 1431 into the solution 
formulas (1311 - (1401) to calculate Uqi an d «02- Inserting these values in (jJTJ) and applying 
the triangle inequality proves the following lemma. 

Lemma 5. Let s = ±i(,o + s' , < \s'\ <C 1, where £o = VPM^ and y?(±z£ ) — with 
< £ < 1- Also assume Re(s') = r] > 0. Then, the solution of (|33|) -( 134|) satisfies the 
boundary estimate 



\u(0)\ < 
\v(0)\ < 



K 

V 
K 



2/i + A 

VP 

2/i + A 

VP 



+ VPI^I 

l^ll + VPI&I 



(45) 
(46) 



where the constant K > is independent of fi and A. TTie solution is as smooth as the 
boundary data and is therefore boundary stable. The solution operator has a simple pole 
at s = ±i£ and as a consequence, the solution in physical space grows linearly in time. 
The growth rate is proportional to \gi\/y/Ji as // — > 0. 
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We shall now discuss the case s — > in more detail. We assume \u\ > u>q > 0, which 
implies s — > 0. Note that the eigenvectors in the general solution (|35|) become linearly 
dependent in the limit, because both K\ = \u\ and k 2 = \oj\ for 5 = 0. We therefore 
assume Re s = fj > 0, and study the the solution in the limit \s\ — > 0. 

Because \s\ <C 1, we can simplify ( 1371) to 

We eliminate «02 using (138]) and obtain 

1 + C ( 1 + TT^W ) - I 1 + - ] I «01 



2 V 2/i + A 




2/i|w| V 2 2/i + XJ 2lu \ 2 
For small |s| 2 we obtain to first approximation 



^ (X + 2fi) 2 g 1 i(\ + 2n)g 2 



s Uo1 = Mi n i nT7 S — • ( 48 ) 



Relation ( |38l) can be written 



S 2 „ 2 #2 

W oi + u 02 = --n 01 - — 



(A + 2/i) 2 ^ | t//ff 2 ^ 



2(A + /x)/x|w| 2(A + /i)w' 

The solution on the boundary is given by (j4T|) . The first component satisfies ix(0) = 
Urn + M02, and ( H9|) shows that u(0) is bounded independently of s. The expression for the 
second component can be simplified for \s\ <C 1. We have to leading order 

-/^\ «M I \ 5 2 \ A z'|u;| / 1 us 2 \ 

i|u;| A z|u;| s 2 f A /i 



(w i + M02) — «oi - : , ^02 • (50) 

w w z \ A + 2/i y 

Therefore, also v(0) is bounded independently of s. The factor in the denominator of 
the right hand side of ( H9l) gives the desired result that our problem is strongly boundary 
stable at s = 0. 



15 



4 Influence of truncation errors on the generalized 
eigenvalues 

Consider the homogeneous differential equations (CQ) with boundary conditions ([3]). Let 

gi = OL\h 2 u xxx + a 2 h 2 v yyy , g 2 = fiih 2 v xxx + (3 2 h 2 u yyy , (51) 

denote the principal part of the truncation error in a second order accurate method with 
grid size h. We can think of boundary conditions ([3]) with boundary data ( 151]) as modified 
homogeneous boundary conditions. Again, we solve the problem using the technique of 
Section [3731 in terms of the simple wave ansatz ( 1351) . By section l3~3| the modified boundary 
conditions become 

(1 - 7 2 )kiK 2 m i + («2 _ 7 2 ^ 2 )mo2 + K 2 gi = 0, 
(tu 2 + K 2 2 )u m + 2u 2 u 2 + iwgz = 0, 

where 



u 4 



gx = -aih (fi^oi + ^2^02) - ot 2 h u KxU i H u 02 

V K 2 

g 2 = i{i x h 2 f^Woi + nluu 02 ^j - i(3 2 h 2 u 3 (u 01 + u m ). 

For small fi, the main effect comes from g%. For simplicity, we therefore assume that 
g 2 = and obtain the equations 

§1 = -h 2 (aiK\ + a 2 u 2 Ki)u i - h 2 ( ai/Co + a 2 — ) u 02 , §2 = 0. (52) 

V K 2j 

Introducing the scaled variable s and the formulas for Kj according to (136]) gives us 
relations (!37|) - (l38|) with g 2 = 0. By using the homogeneous equation ( 138|) . we can eliminate 
u 02 from (137|) and (152]) . resulting in the solution formula ( 139]) with g 2 = 0, and 

gi = -h 2 (aiK\ + a 2 uj 2 Ki)u m + h 2 ^1 + + "2^-^ Uoi- 



Hence, the solution formula (1391) defines a perturbed eigenvalue problem that can be 
written in the form 

(f(s)u i = d(s 2 )u m . 

Since 

K 2 PfJ, 



\u\ \ ' A + 2/i' 
16 



we have 



~ 2 x (A + 2/i)/i 2 K 2 ( 3 2 / S 2 \ / 3 W 4 



) = — ?r-rn — r-r + «i - 1 + — a x td + a 2 — 

2fJL\U\ \U)\ \ \ 2 J \ K 2 

We assume now that that A//i ^> 1. For the unperturbed problem, the properties of 
the generalized eigenvalues are given in Table [TJ 

s 2 ~ -0.9, ^4 ~ 0.3, 



a; a; 



Therefore, 



6 (si) ~ fa 1T ^M 3 + a 2 w 2 ^|u;| - 0.55 fa 1T ^-|w| 3 + a 2 w 2 M— 11 

2/x|a;| |w| \ \u>\ 6 \u>\ \ \u}\ 6 n 2 J J 

\h 2 u 2 

~ — (0.027a! + 0.3a 2 - 0.55(ati + a 2 )). 

We want to estimate how sensitive the generalized eigenvalues s = ±z£ are to trun- 
cation error perturbations. We perturb (p(s) around Sq. For A//i > 1, Table [Hand ( I30p 
gives 

^^^M^4-2fl + ^V-0.67 
\u)\ | a; | s w \ 2 J 

„ , 0.67 
p'(s ) ~s -^y~±2.12 ? . 

The Taylor expansion ( 1421) gives for small /icj, 

(5 - SoV(so) = 0(3o)- 

We get 

flfs2^ \h 2 uj 2 i 
s-s ~ -y*L ~ T ^_( .027 ai + 0.3a 2 - 0.55( ai + a 2 )) — -. (53) 
</? '(so) 2/i 2.12 

We now make some observations. Because 9(sq) is real, the generalized eigenvalue is 
perturbed along the imaginary axis and remains purely imaginary. Hence the perturbed 
problem is well-posed. The value of the perturbed generalized eigenvalue determines the 
phase velocity of surface waves in the numerical solution. To avoid large phase errors, we 
must therefore keep the perturbation of the generalized eigenvalue small. If we accept a 
relative error in the phase speed of size e, where < e < 1, we have to choose the grid 
size h such that 

\h 2 u 2 \a \ 



e « 1. (54) 
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If the computational grid has P grid points per wave length L = 2tt/\u\, we get 

P wP' 1 1 P \ e nj 



Hence, the number of grid points per wave length must be proportional to yXjJi to 
maintain the accuracy as /i/A — > 0. 

For a fourth order accurate method, where the leading order truncation error terms 

are 

^ = ^— + a 2 h — , 

equation (I54I) is replaced by 

\ *,4, ,4| / I 

^1 = e « 1. (55) 



\a n A x 



/i 

The number of grid points per wave length to maintain an e-error in the phase velocity 
now becomes 

Therefore, as /i/A — > 0, the number of grid points per wave length grows much slower for 
the 4th than the 2nd order accurate method. 

For other truncation error perturbations of the boundary conditions, such as a u xxxx 
term in g\, 6(§l) becomes complex. If the truncation error coefficient has the wrong sign, 
the perturbed problem gets eigenvalues with positive real part. From Lemma [1] we know 
that such problems are ill-posed. Furthermore, the factor /i in the denominator of ( 153]) 
shows that the rate of the exponential growth can get arbitrarily large as yu/A — > 0. It is 
therefore very difficult to compensate for such growth with an artificial dissipation term. 



5 Numerical experiments 

For a second order hyperbolic equation, energy conservation ensures that all eigenvalues 
of the spatial operator are either real and negative, or zero. The same property applies 
to the discretized problem. To avoid any spurious growth in the numerical solutions, 
it is therefore important to use a discretization that also satisfies energy conservation. 
Such a discretization was derived for the 3-D elastic wave equation in Nilsson et al [?]. 
In the present work, we use the corresponding discretization for the two-dimensional 
case. This numerical method discretizes the elastic wave equation with a second order 
accurate, energy conserving, finite difference method on a Cartesian grid with constant 
grid sizes in space and time. The second order method was recently generalized to fourth 
order accuracy by Sjogreen and Petersson [8], and we use both the second and forth 
order methods in the following numerical experiments. Note that our finite difference 
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methods are based on solving the elastic wave equation as a second order hyperbolic 
system using summation by parts operators. These methods are fundamentally different 
from the commonly used staggered grid method developed by Vireaux [9], Levander [6], 
and others, which is based on solving the elastic wave equation as a first order hyperbolic 
system. 

5.1 Surface waves 

To study surface waves using real arithmetic, we are interested in the real part of the 
eigenfunction fl32l) corresponding to the generalized eigenvalue 

s = i£ - 

Assuming oj > 0, the real part of ( 1321) can be written as 

fT~t2 I cos ( u(y + c r t)) 

\y/l-$sin(u;(y + c r t)) j 

+ e Wi-^/(HA), ( cos(u(y + c r t)) 



2 



sm 



(u(y + c r t))/ v /l-4 2 /i/(2/x + A) / 



Here, we define the Rayleigh phase velocity by 

c r — io\/J*>- 

To perform reliable numerical simulations, it is of great interest to know the number of 
grid points per wave length, P, that is required to obtain a certain accuracy in a numerical 
solution. If the wave length is L = 2tt/\uj\, we define 

L 

We consider a periodic domain in the ^-direction and choose the computational domain 
to contain exactly one wave length of the solution. In this investigation we shall keep the 
wave length fixed at L — 1, which gives the spatial frequency uo = 2ir. For simplicity, 
we set A = 1 in all numerical experiments. A free surface boundary condition is imposed 
at x = 0. We truncate the computational domain at x = L x where we impose an 
inhomogeneous Dirichlet condition. The boundary data is given by the exact solution 
( 1561) . which is exponentially small along x = L x . For all values of X/fi (see Table [Q, 
«i/|o;| > 0.2955, and we make the influence of the Dirichlet boundary closure small by 
choosing 

L x = 10, e -^V^I^ < e-27T-0.2955.10 ^ 8 5 . 1Q -9 
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2nd order method 4th order method 




time time 

Figure 3: Max error as function of time for a Rayleigh surface wave with A = 1 and 
H = 0.01. Results from the second and fourth order methods are shown on the left and 
right, respectively. The different colors corresponds to different number of grid points per 
wave length. Note that the grids are coarser for the fourth order computations. 



In our first experiment, we take [l = 0.01. The numerical solution is evolved from 
initial data given by (I56I) at time t = and t = —5 t , where the time step satisfies the 
Courant condition (recall that we have scaled time to give unit density) 



K c === , K c 




second order method, 
\/A + 3/i' | 1.3, fourth order method. 



In Figure [3] we show the max norm of the error in the numerical solution as function of 
time for t < 20. Since the wave length in the ^-direction is one, the number of grid points 
per wave length satisfies P — N y — 1. Results for the second order accurate method are 
shown on the left, illustrating the expected convergence rate as the grid is refined. Note 
that at least 100 grid points per wave length (green line) are needed to obtain a numerical 
solution to within about 5% of the exact solution. On the right side of the same figure, 
we show results for the fourth order method. Here the error decreases by a factor of 16 
when the number of grid points is doubled. In this case, only 20 grid points per wave 
length are needed to make the error less than about 5% of the exact solution. 

In our next experiment, we study how the accuracy depends on /i when the second 
order method is used for propagating the Rayleigh wave (1561) . The period of the wave is 

T=*L = ! = J^. (57) 

UC r C r ^Oa/A* 
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Case 
H = 0.1 
T = 3.330 



P = 1/h 
25 
50 
100 



\u ( , 



o(t = T) 
5.26 • 10~ 2 
1.48 • 10~ 2 
3.85 • 10~ 3 



,(t = 10T) 



2.99 • 10" 1 
9.26 • 10~ 2 
2.44 • 10~ 2 



n = 0.01 
T = 10.474 



25 
50 
100 
200 



2.45 • 10" 



1.07-10 



-i 



3.32 • 10~ 2 
8.86 • 10~ 3 



7.64 • 10" 1 
5.59 • lO" 1 
2.06 • 10" 1 
5.73 • 10~ 2 



n = 0.001 
T = 33.104 



100 
200 
400 



1.96 • 10 



-i 



7.40 • 10~ 2 
2.13 • 10~ 2 



7.69 • lO" 1 
4.25 • lO" 1 
1.36 • 1Q- 1 



Table 2: The max norm of the error in the numerical evolution of the Rayleigh surface 
wave, after one and ten periods. Note how the number of grid points per wave length, 
P = 1/h, must be drastically increased to maintain the accuracy as /J becomes smaller. 



In Table |5] we show the max norm of the error after one and ten periods. Note that the 
period gets longer, i.e., the surface wave propagates slower as \i — > 0. The case \i = 0.1 
shows close to second order convergence, both at time t = T and t = 10 T. The error 
levels are reasonable for a second order method, but increase with time because the error is 
dominated by phase errors, i.e., the numerical solution propagates with a slightly different 
phase velocity compared to the exact solution. The error gets larger for \i = 0.01, and a 
finer grid must used to obtain comparable error levels. For /x = 0.001, the grid must be 
refined further to obtain reasonable error levels, and the cases P = 25 and P = 50 are 
inadequate. A visual inspection shows that after 10 periods, the numerical solution with 
P = 50 is more than 180° out of phase with the exact solution (experiment not shown to 
save space). We only observe close to second order convergence when the grid is refined 
from 200 to 400 grid points per wave length. 

Note that for /j = 0.001, the grid with 200 grid points per wave length gives of the 
order 10 percent accuracy after one period (HugHoo ~ 0.545). This grid is about 10 times 
finer than what is normally required to get that accuracy with a second order method [1] . 
In the a;-direction, the gradient of the exact solution is the largest along x — 0, and 
l^xl = \u y \ for all \x. In the limit \i — > 0, it is straight forward to show \u x \ = \v y \. Hence, 
the gradient of the exact solution is of the same order in both directions, and conclude 
that solution is extremely well resolved on the grid. Furthermore, the phase velocity of 
the surface wave becomes slower and slower as /i — > 0, while the time step is governed 
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by y/X + 3fi, which tends to \/A = 1. Hence, the temporal resolution of the surface wave 
only improves as — > 0. 

The analysis of the phase velocity in §|4] shows that truncation errors in a second order 
accurate method perturb the generalized eigenvalue according to 

i = £o + e, e = • (58) 

/i 

The perturbed generalized eigenvalue corresponds to a perturbed phase velocity c' r = y/Ji^. 
Assuming that phase errors dominate the numerical errors, the amplitude of the error 
follows by 

e(t) = u(c' r — c r )t = ojy/Jlet. 

The period of the surface wave follows from (I5T1) . so e(T) = C\eT, C\ = const. For a 
computational grid with grid size h = 1/P, fl58l) gives 

6 = PV 2 = 

Hence, to maintain a constant error level in the numerical solution after a fixed number 
of periods, we must choose Py/JI = const., if A is constant. This assertion is tested by 
the numerical experiment shown on the left side of Figure HI Here we show the max 
error as function of time scaled by the period of the solution. The first case (red curve) 
corresponds to /i = 0.1, with period T = 3.33 and resolution P = 40 grid points per 
wave length. Notice how closely this error curve follows the case /i = 10~ 3 , with period 
T = 33.104 and a grid with 400 grid point per wave length. We conclude that the second 
order method needs a prohibitively fine computational grid to accurately calculate surface 
waves for small values of fi. 

We repeat the above experiment with a fourth order accurate method. The results 
are shown on the right side of Figure HI In this case we obtain similar error levels using 
a significantly coarser grid. For \i — 0.1 and fi = 0.001, we use P = 12 and P = 38, 
respectively. For the fourth order method, the perturbation of the generalized eigenvalue 
is given by fl55|) . Using the same argument as for the second order method, we must 
choose P/i 1 / 4 = const, to obtain a constant error level in the numerical solution after a 
fixed number of periods. This scaling is approximately preserved in these calculations, 
since 

'6.748, P = 12, (j, = 0.1, 
6.757, P = 38, fi = 0.001. 



P/i 



1/4 



We conclude that the fourth order method is much better suited for simulations when /i 
is small. Compared to the second order method, the fourth order method needs a smaller 
number of grid points per wave length, and the required resolution grows much slower as 
H ->• 0. 
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Figure 4: Max error in the numerical evolution of the Rayleigh surface wave, as function 
of time scaled by the period, T = l/(£o-y/p)- For the second order method (left), the case 
\l = 0.1 with P = 40 is shown in red and fi = 0.001 with P = 400 is shown in blue. For 
the fourth order method (right), the case /i = 0.1 with P = 12 (green) and ji = 0.001 
with P = 38 (black) give comparable error levels. 

To indicate how much more efficient the fourth order method is in practice, we give 
some execution times obtained on a MacBook Pro laptop computer. The above numerical 
experiments for /x = 0.001 required 20, 604 seconds (~ 5 hours, 43 minutes) for the second 
order method with P = 400. Similar accuracy was obtained with the fourth order method 
using P = 38, but this calculation only took 60 seconds. Hence, for this problem the fourth 
order method was 343 times faster than the second order method. 

5.2 Mode to mode conversion 

Consider a compressional wave of unit amplitude traveling in the negative x-direction in 
a homogeneous material, with displacement 




e m+kx+uy) ^ £ = COS0>O, w = s in0, £>0. 
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If this wave encounters a free surface boundary at x = 0, it will be reflected and split into 
two waves that both travel in the positive x-direction, 



u (oni) =u (P )+u (5) ) 

u< p > =R p ( ~ k \ e ^- fc ^), 

U (S) = gj f ] f /( t /-.A, + ... v , n ;> 

Va 2 k 2 + u 2 y-aife 

The reflected waves correspond to a compressional and a shear wave, since the curl of 
u( p ) and the divergence of are zero. In order for vft 71 ' and u* "*) to satisfy the elastic 
wave equation (CQ) with F\ = F 2 = 0, the frequency and wave numbers must satisfy the 
elementary relations 

£ 2 = (A + 2/i)(A; 2 + w 2 ) = A + 2/x, £ 2 = /i(« 2 A; 2 + oo 2 ). (59) 

We select the signs of £ and a such that ir m ) and u^*-* travel in the negative and positive 
x-direction, respectively. The amplitudes of the reflected waves, R p and R s , are functions 
of A, [A, and the angle of the incident wave, <p. The amplitudes Rp and R s are uniquely 
determined by the free surface boundary conditions ([3]) (with g\ — g 2 — 0). For a more 
detailed discussion, we refer to Achenbach [T], § 5.6. 
As a consequence of the relation (1391 . 



a 2 = l + - A + /i 



cos 2 . 



Hence, when /i < A, the reflected S-wave will propagate almost parallel to the x-direction 
because a 2 ^> 1, see Figure [5j The wave lengths of the compressional and shear waves 
are given by 

2tt / // 

L„ = — : = 27T, L s = 2lX a - . 

P Vk 2 + u 2 y \ + 2fi 

Note that the wave length of the compressional wave is fixed, while L s becomes small as 

To include two wave lengths of u^ n ) in the computational domain, we take L y = 
An/ sin</> and L x = An/ cos<p. As before, we impose periodic boundary conditions in 
the ^/-direction, a Dirichlet boundary condition at x = L x and a free surface condition 
at x = 0. By construction, the function u < - m - ) + u* "*) is L^-periodic in the y-direction, 
satisfies the elastic wave equation in the interior, and the free surface condition at x = 0. 
In principle, we could compute a numerical approximation of u*" 1 ) + u* "*) by adding a 
suitable forcing function to the Dirichlet boundary condition at x = L x . However, we 
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2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 



Figure 5: f-component of the outgoing shear wave as function of (x,y) at t = 0. The 
angle of the incoming P-wave is = tt/4. The frames correspond to \x = 1.0 (left), \x = 0.1 
(middle), and fi = 0.01 (right). 



instead choose to only compute the outgoing S-wave, u^ s \ For this reason, we impose the 
inhomogeneous Dirichlet boundary condition 

u(L x ,y,t) = u (s) (L x ,y,t), 

and take the forcing functions in the normal stress boundary conditions ([3]) to be 

Si = -W m) +4 P) )-7 2 {v$*+v<F), 

92 = - (4 m) + 4 P) + + ^ ] ) ■ 

We use the exact solution as initial conditions for the numerical solution. 

To accurately solve this problem numerically, it is necessary to resolve the short shear 
waves on the computational grid. For this problem, we define the resolution in terms of 
the number of grid points per shear wave length, 

h h VA + 2/i' 

We evaluate the error in the numerical solution as function of time for two materials. 
The first material has (A = 1, \x = 0.1) and the second has (A = 1, \x = 0.01). As a 
consequence, the period of the wave is slightly different for the two cases 

t _2tt 2tt ^ f 5.74, ft = 0.1, 

" T ~ VA + 2/x ~ [6.22, fi = 0.01. 

In Figure[6]we show the error as function of normalized time, t/T, for the two materials, 
using the fourth order accurate method. Note that the error levels are comparable for 
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Figure 6: Results for computing the outgoing shear wave with different resolution, char- 
acterized by the number of grid points per wave length, P s . The relative error in max 
norm is shown as function of time scaled by the period of the wave. Two cases are shown, 
(A = 1, fi = 0.1) and (A = 1, n = 0.01). 

the same number of grid points per wave length, and converge to zero as 0{P~ i ) as 
the grid is refined. Thus the mode to mode conversion problem does not suffer from 
the same extreme resolution requirements as the surface wave problem in the previous 
section. Because we have scaled the problem such that the P- waves have wave length 2ir, 
the S-waves get a wave length of the order 27r^//J. Hence, to keep the number of grid 
points per S-wave length constant for different materials, we have to choose the grid size 
according to 

h = 2n ^ 
+ P s ' 

Compared to the material with /i = 0.1, the grid size must therefore be taken about a 
factor of vlo smaller for the case jj, = 0.01, to obtain the same number of grid points 
per wave length. This scaling is independent of the order of accuracy in the numerical 
method. 

No surface waves can be triggered by a propagating P-wave because the relation ( 159]) 
shows that £ 2 /(/^u; 2 ) > 1. However, evanescent modes due to an interior forcing function 
could trigger both S-waves and surface waves. Since the surface waves are only slightly 
slower than the S-waves, their wave length is of the same order as the length of an S-wave 
of the same frequency. If the problem is scaled such that the P-wave length is constant, 
both the S-wave and the surface waves would therefore have wave lengths of the order 




mu=0.1 , P=20 
mu=0.1 , P=40 
mu=0.01, P=10 
mu=0.01 , P=20 
mu=0.01 , P=4Q 
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JJL. Based on the results of Section I5.1[ a second order accurate method would need a 
grid size of the order h ~ /j, to maintain a constant accuracy in the numerical solution as 
fj, — > 0. For a fourth order method, it would suffice to use h ~ /i 3//4 . 

6 Conclusions 

We have developed a normal mode analysis for the half-plane problem of the elastic wave 
equation subject to a free surface boundary condition. Our analysis allows the solution 
to be estimated in terms of the boundary data, showing that the solution is as smooth 
as the boundary forcing. Hence, using the terminology of [5], the problem is boundary 
stable. The dependence on the material properties is transparent in our estimates. Using 
a modified equation approach, the normal mode technique was extended to analyze the 
influence of truncation errors in a finite difference approximation. Our analysis explains 
why the number of grid points per wave length must be so large when calculating surface 
waves in materials with /i/A 1. To obtain a fixed error in the phase velocity of Rayleigh 
surface waves, our analysis predicts that the grid size must be proportional to /i 1 / 2 for a 
second order method, when A = const. For a fourth order method, the analysis shows that 
it suffices to use h ~ /i 1//4 . These scalings have been confirmed by numerical experiments. 

It is theoretically possible to derive stable finite difference schemes that give higher 
than fourth order accuracy. These methods use wider stencils that are more expensive 
to evaluate, but for the surface wave problem, it would suffice to use a grid size of the 
order h ~ yU 1,/p , where p is the order of accuracy. For sufficiently small values of [i these 
methods should be more efficient as the order of accuracy increases. However, numerical 
experiments must be performed to evaluate how small fi must actually be to compensate 
for the higher computational complexity of these very high order accurate methods. 
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A Miscellaneous lemmata 

Lemma 6. Let u be a real number and let s = rj + it; be a complex number where rj > 0. 
Consider the relation 

k = y/oo 2 + s 2 , 

where the branch cut in the square root is defined by 

— 7r < arg(a + i(3) < n, arg \J a + ij3 = - arg(a + i(3). (60) 
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Then, 

Re(K)>rj. (61) 

Proof. Since r\ > 0, we can write 

s = 77(l + i£'), /t = r/A/o; /2 + (l+^e') 2 ) £' = ~, w' = -. (62) 

v V rj 

Define real numbers a and b such that 



a + i& = V& 2 + (1 + i£') 2 , a ^°- ( 63 ) 
Squaring relation fl63|) and identifying the real and imaginary parts give 

oh = 

2 l2 1 f/2 . /2 

a — 6=1 — 4 + u . 
The first relation gives b = £'/ a, which inserted into the second relation results in 



a 2 -^ = l-£' 2 + u/ 2 . (64) 



Note that the left hand side is a monotonically increasing function of a 2 . When u>' = 0, 
equation ( )64l) is solved by a 2 = 1. The right hand side of (|64l) is a monotonically increasing 
function of u' 2 . Therefore a 2 > 1 for a/ 2 > 0. We conclude that the unique solution of 
(EH) satisfies 



a 2 > 1. 



Because a must be non- negative, we have a > 1. Relations (|62|) and ( 163]) give 



Re («) = rj Re y u/ 2 + (1 + z£') 2 = 77 a > rj. 

□ 

Corollary 1. Let u be a real number, s = rj + i£ be a complex number with rj > 0, and 
let < 7 < 00 be a constant. Then there is another constant < 5 < 00 such that 



Re\juj 2 + ^ > <fy, r/ = Re(s)>0. (65) 



Proof. Let s' = s/7. Lemma |6] proves that 



ReJu 2 + ^ = Re v 7 ^ 2 + s' 2 > Re (s') = -Re (s). 
V 7 7 

Hence, 5 = I/7 > and the corollary follows. □ 
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B The case uj -> 

We now extend the boundary estimates in Section 13.31 to the case u> — > when s = sq, 
| so| > is fixed. In this limit, 



U~ 2 I 



For \s\ 3> 1, we can simplify ( 1371) according to 



s 2 



s 2 



OO. 



2fi\u\ \ X + 2/j, 



sgi- 



In a similar way, ( 1381) becomes 



2 



^92- 



2ujs 2 

Solving the latter equation for w i and inserting into (166]) gives 

(A + 2/i)s^i 



s 2 



For large \s\, we have to leading order, 



2fi\u\ 



V ,192 



A + 2/i 2u V A + 2/x 



U 2 = -i 



A + 2/i 



A + 2/i cus 2 



Note that s = sJJI\oj\, and 

1 



Hence, 



S\0J\ 

and (jSTJ give 



s 



s 2 M 



.s 2 



«02 



^ 92 

2s 2 



uj -»■ 0. 



The solution on the boundary follows from (f4T|) and gives directly 

vrrt - I a VA + 2/i^! 

■U(O) = M 1 + M02 = h U{0J). 



For large |s|, we can simplify the expression for t>(0), 

fi(0) = 



i|a>|s«oi «wuo2 /A + 2/4 \f\ig2 



UJ 



\UJ\S 



2s 



+ 0{uj) 



(66) 



(67) 



(68) 



The expressions for u(0) and -0(0) show that the solution is well behaved in the limit 
oj ->■ 0. 
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